{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5d904dee",
   "metadata": {},
   "source": [
    "# Example 6: Solving Partial Differential Equation (PDE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d568912",
   "metadata": {},
   "source": [
    "### We aim to solve a 2D poisson equation $\\nabla^2 f(x,y) = -2\\pi^2{\\rm sin}(\\pi x){\\rm sin}(\\pi y)$, with boundary condition $f(-1,y)=f(1,y)=f(x,-1)=f(x,1)=0$. The ground truth solution is $f(x,y)={\\rm sin}(\\pi x){\\rm sin}(\\pi y)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0e2bc449",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pde loss: 5.26e+00 | bc loss: 7.96e-02 | l2: 2.62e-02 : 100%|█| 20/20 [00:19<00:\n"
     ]
    }
   ],
   "source": [
    "from kan import KAN, LBFGS\n",
    "import torch\n",
    "import matplotlib.pyplot as plt\n",
    "from torch import autograd\n",
    "from tqdm import tqdm\n",
    "\n",
    "dim = 2\n",
    "np_i = 21 # number of interior points (along each dimension)\n",
    "np_b = 21 # number of boundary points (along each dimension)\n",
    "ranges = [-1, 1]\n",
    "\n",
    "model = KAN(width=[2,2,1], grid=5, k=3, grid_eps=1.0, noise_scale_base=0.25)\n",
    "\n",
    "def batch_jacobian(func, x, create_graph=False):\n",
    "    # x in shape (Batch, Length)\n",
    "    def _func_sum(x):\n",
    "        return func(x).sum(dim=0)\n",
    "    return autograd.functional.jacobian(_func_sum, x, create_graph=create_graph).permute(1,0,2)\n",
    "\n",
    "# define solution\n",
    "sol_fun = lambda x: torch.sin(torch.pi*x[:,[0]])*torch.sin(torch.pi*x[:,[1]])\n",
    "source_fun = lambda x: -2*torch.pi**2 * torch.sin(torch.pi*x[:,[0]])*torch.sin(torch.pi*x[:,[1]])\n",
    "\n",
    "# interior\n",
    "sampling_mode = 'random' # 'radnom' or 'mesh'\n",
    "\n",
    "x_mesh = torch.linspace(ranges[0],ranges[1],steps=np_i)\n",
    "y_mesh = torch.linspace(ranges[0],ranges[1],steps=np_i)\n",
    "X, Y = torch.meshgrid(x_mesh, y_mesh, indexing=\"ij\")\n",
    "if sampling_mode == 'mesh':\n",
    "    #mesh\n",
    "    x_i = torch.stack([X.reshape(-1,), Y.reshape(-1,)]).permute(1,0)\n",
    "else:\n",
    "    #random\n",
    "    x_i = torch.rand((np_i**2,2))*2-1\n",
    "\n",
    "# boundary, 4 sides\n",
    "helper = lambda X, Y: torch.stack([X.reshape(-1,), Y.reshape(-1,)]).permute(1,0)\n",
    "xb1 = helper(X[0], Y[0])\n",
    "xb2 = helper(X[-1], Y[0])\n",
    "xb3 = helper(X[:,0], Y[:,0])\n",
    "xb4 = helper(X[:,0], Y[:,-1])\n",
    "x_b = torch.cat([xb1, xb2, xb3, xb4], dim=0)\n",
    "\n",
    "steps = 20\n",
    "alpha = 0.1\n",
    "log = 1\n",
    "\n",
    "def train():\n",
    "    optimizer = LBFGS(model.parameters(), lr=1, history_size=10, line_search_fn=\"strong_wolfe\", tolerance_grad=1e-32, tolerance_change=1e-32, tolerance_ys=1e-32)\n",
    "\n",
    "    pbar = tqdm(range(steps), desc='description')\n",
    "\n",
    "    for _ in pbar:\n",
    "        def closure():\n",
    "            global pde_loss, bc_loss\n",
    "            optimizer.zero_grad()\n",
    "            # interior loss\n",
    "            sol = sol_fun(x_i)\n",
    "            sol_D1_fun = lambda x: batch_jacobian(model, x, create_graph=True)[:,0,:]\n",
    "            sol_D1 = sol_D1_fun(x_i)\n",
    "            sol_D2 = batch_jacobian(sol_D1_fun, x_i, create_graph=True)[:,:,:]\n",
    "            lap = torch.sum(torch.diagonal(sol_D2, dim1=1, dim2=2), dim=1, keepdim=True)\n",
    "            source = source_fun(x_i)\n",
    "            pde_loss = torch.mean((lap - source)**2)\n",
    "\n",
    "            # boundary loss\n",
    "            bc_true = sol_fun(x_b)\n",
    "            bc_pred = model(x_b)\n",
    "            bc_loss = torch.mean((bc_pred-bc_true)**2)\n",
    "\n",
    "            loss = alpha * pde_loss + bc_loss\n",
    "            loss.backward()\n",
    "            return loss\n",
    "\n",
    "        if _ % 5 == 0 and _ < 50:\n",
    "            model.update_grid_from_samples(x_i)\n",
    "\n",
    "        optimizer.step(closure)\n",
    "        sol = sol_fun(x_i)\n",
    "        loss = alpha * pde_loss + bc_loss\n",
    "        l2 = torch.mean((model(x_i) - sol)**2)\n",
    "\n",
    "        if _ % log == 0:\n",
    "            pbar.set_description(\"pde loss: %.2e | bc loss: %.2e | l2: %.2e \" % (pde_loss.cpu().detach().numpy(), bc_loss.cpu().detach().numpy(), l2.detach().numpy()))\n",
    "\n",
    "train()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2246bab",
   "metadata": {},
   "source": [
    "### Plot the trained KAN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "02e2a0ba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABchElEQVR4nO3dd1yVZf8H8M99DnsoUxRBGaICoqi40nJrucssFXA1HaFmjtRylLM098rKBMydlpia28rBdLBEwAGi7L3OOff1+4Pnvn9QDsaBs77v1+v551E4l51znc99re/FMcYYCCGEECWSqLoBhBBCtA+FCyGEEKWjcCGEEKJ0FC6EEEKUjsKFEEKI0lG4EEIIUToKF0IIIUpH4UIIIUTpKFwIIYQoHYULIYQQpaNwIYQQonQULoQQQpSOwoUQQojSUbgQQghROgoXQgghSqen6gYQogkYY8jKykJhYSHMzMxgbW0NjuNU3SxC1BaNXAh5gdzcXGzcuBFubm6wtbWFs7MzbG1t4ebmho0bNyI3N1fVTSRELXF0EyUhz3b69GmMHj0axcXFACpGLwJh1GJiYoIjR45g8ODBKmkjIeqKwoWQZzh9+jSGDh0Kxhh4nn/u35NIJOA4DiEhIRQwhFRC4ULIv+Tm5sLBwQElJSUvDBaBRCKBsbExUlJSYGFhUf8NJEQD0JoLIf/y888/o7i4uFrBAgA8z6O4uBh79+6t55YRojlo5EJIJYwxuLm5ISkpCTXpGhzHwcXFBQkJCbSLjBBQuBBSRWZmJmxtbev089bW1kpsESGaiabFCKmksLCwTj9fUFCgpJYQotkoXAipJCsrq04/b25urqSWEKLZKFyIzsvJycGuXbvQu3dvdO3atdZrJhKJBF988QWuXbtWo/UaQrQRhQvRSWVlZTh27BhGjx4Ne3t7TJ8+HSYmJggMDMTq1atrFTB9+vTB8ePH0aNHD7Ru3RrLli3DvXv36qH1hKg/WtAnOoMxhr///hvBwcE4dOgQcnJy0LFjR/j5+WHs2LFo2rQpgLqdczE3N8fFixcRGBiII0eOoLCwED169ICfnx/effddWuwnOoPChWi9+Ph4BAcHY9++fUhOToajoyN8fX0xfvx4eHp6PvNnanpC/+TJkxg0aFCVPysuLsbx48cRFBSE06dPQyKR4I033oC/vz+GDRsGIyMjpf47CVEnFC5EK6Wnp+PgwYMICgpCaGgoGjVqhLfffht+fn549dVXIZG8fEa4urXFjh49+p9g+benT5/iwIEDCAwMRFhYGBo3bowxY8bUqD2EaBIKF6I1iouL8fvvv4sjBY7j8MYbb8DPzw9Dhw6FsbFxjX9nbm4u9u7di02bNiExMVH8/11dXREQEICJEyeicePGNfqdcXFxCAoKQlBQEB48eIAWLVrA19cX/v7+cHd3r3EbCVFHFC5EoykUCly+fBmBgYE4evQoCgoK0L17d/j5+WHMmDGwsbFRyuswxpCdnY2CggKYm5vDysqqzifxeZ7HP//8g8DAQBw8eBC5ubno1KkT/P39q6wBEaKJKFyIRrp9+7a4jpKamgpXV1f4+vrC19cXrVq1UnXzaqysrAwhISEICgrCiRMnoFAoMGjQIPj5+WHUqFEwNTVVdRMJqREKF6IxHj9+jF9++QXBwcG4efMmrKys8O6778LX1xfdu3fXmppe2dnZOHToEAIDA/H333/D1NQUb731Fvz9/dGvXz9IpVJVN5GQl6JwIWqtoKAAx44dQ2BgIM6fPw8DAwMMGzYM/v7+GDx4MAwMDFTdxHqVlJSE4OBgBAYGIiEhAc2aNcP48ePh5+eHDh06aE2gEu1D4ULUjlwux9mzZxEUFIRjx46hpKQEvXv3hq+vL0aPHq2Td6YwxhAaGoqgoCD88ssvyMzMRLt27eDn5wdfX184ODiouomEVEHhQtQCYwzh4eEIDg7G/v37kZ6eDnd3d/j5+WHcuHFo2bKlqpuoNmQyGc6cOYPAwEAcP34cZWVl6NOnD/z9/TF69Gg0atRI1U0khMKFqNb9+/exb98+BAcHIy4uDnZ2dhg7diz8/PzQsWNHmvZ5ifz8fBw5cgRBQUG4cOECDA0NMWLECHHaUF9fX9VNJDqKwoU0uJycHPEL8cqVKzAxMcGoUaPg5+eH/v37Q09PT9VN1EiPHj3CL7/8gsDAQNy5cwc2NjZiUNelICchtUHhQhpEeXk5/vjjD3GrrVwuR79+/eDv749Ro0bBzMxM1U3UGowx3Lp1C4GBgdi3bx/S0tLg5uYGPz8/+Pn5wcXFRdVNJDqAwoXUG8YYrl69iqCgIBw6dAjZ2dnw9vaGr68vxo4dC3t7e1U3UespFAqcP38eQUFBOHLkCIqKitCzZ0/4+fnhnXfegZWVlaqbSLQUhQtRuoSEBAQHByM4OBhJSUlwcHDA+PHj4evri3bt2qm6eTqrqKgIx48fR2BgIM6cOQOpVIqhQ4fC398fQ4cOhaGhoaqbSLQIhQtRioyMDBw8eBDBwcG4fv06zM3NMXr0aPj7++O1116jwoxq5smTJ9i/fz+CgoIQHh4OCwsLvPPOO/Dz80PPnj3p/SJ1RuFCaq2kpAQnTpxAYGAgTp8+DQAYPHgw/Pz8MHz48FoViiQNLzY2Viyk+fDhQzg5OcHX1xd+fn5o27atqptHNBSFC6kRnudx+fJlcQ4/Pz8f3bp1g6+vL9555x3Y2tqquomklniex19//YXAwEAcOnQIeXl58PHxEc8aNWnSRNVNJBqEwoVUS3R0NIKCgrBv3z6kpKTAxcVFvHCrdevWqm4eUbLS0lKEhIQgMDAQJ0+eBM/zGDRoEPz9/TFy5EiYmJiouolEzVG4kOdKS0sT5+WjoqJgaWkpzsv36NGDzk3oiKysLBw8eBCBgYG4evUqzMzMMHr0aPj5+aFv375USJM8E4ULqaKwsBDHjh1DUFAQzp07Bz09PQwbNgx+fn54/fXXaUeRjktMTBQLad67dw/29vYYP348/P390b59e1U3j6gRChcCuVyO8+fPIzAwEMeOHUNxcTFeffVV+Pn5YfTo0bC0tFR1E4maYYzhxo0bCAwMxP79+5GVlQUvLy/4+/tj/PjxaN68uaqbSFSMwkVHMcYQGRkpFop88uQJ2rRpAz8/P4wfPx5OTk6qbiLREOXl5Th9+jSCgoJw/PhxlJeXo1+/fuLDibm5uaqbSFSAwkXHPHjwQLxwKyYmBk2aNMHYsWPh6+uLzp070zoKqZO8vDwcOXIEgYGBuHjxIoyNjTFy5Ej4+/tj4MCBVEhTh1C46IDc3FwcPXoUQUFBuHTpUpUO379/f+rwpF48fPgQ+/btQ2BgIGJiYmBra4tx48bBz88PPj4+9CCj5ShctFTlqYrff/+9ylTFm2++SVMVpMEwxhAVFSVuZa88Bevr6wtnZ2dVN5HUAwoXLcIYw/Xr1xEUFISDBw8iKysL7du3h5+fH8aOHUuLrETlKm8eOXr0KIqLi9GrVy/4+/tjzJgxtHlEi1C4aIF79+6JF25V3h7q6+tL20OJ2qq87f3PP/+ssu19yJAhtO1dw1G4aKjMzEwcOnQIQUFBuHbtWpWDbb1796aDbUSjCAd2AwMDERkZWeXAbs+ePWl9RgNRuGiQ0tJSnDhxAkFBQfjjjz/AGBNLcgwfPpxKchCtIJQaCg4OxqNHj+Ds7AxfX1/4+/tTqSENQuGiARhjmDFjBn755Rfk5eWhS5cu8PX1xbvvvkvFBInWqlwk9dChQ8jPz0fXrl2xbt069OrVS9XNIy9B4aIheJ4HAHAcR1MEROcIX1OMMeoDGoLChRBCiNLpqboBmkbXs5ieGAn1AeoD1UHhUgshISGIjIzUqQ8ZYwwDBgxAjx49VN0UogZCQkIQERGhc31g4MCB1AeqicKlFk6dOoWmTZvC3d292j8jrJlo6t3kFy9exLVr16hjEQAVfcDOzg4eHh4N8noKhQISiUSlYXbx4kVcvXqV+kA1UbjUAsdxeO211/Daa6+99O8qFAqcOnUKe/fuhVQqxcSJEzFgwACNO4dSUFCAnJwcVTeDqAmO49C7d2/07NkTMpkMMpkMZmZmSv/yLy4uxg8//IAzZ86gQ4cOCAgIgK2trUpCJj8/n/pADVC41CPGGA4cOICpU6eisLAQAHDy5Els3rwZfn5+OjWlQLSPTCZDQEAArl27hsaNG+PXX39F48aNlfb7FQoF1q5di5UrV0Iul+PkyZO4c+cOAgMDqTaeBtDMORoNcf/+fSxevBiFhYWQSCSQSCQoKCjA559/jpiYmAZdGOV5HvHx8Th16hSSk5N1flGW1B3Hcbh37x4iIyMRExMjPkApA2MMd+7cwZYtWyCXy8X/LyQkBMHBwfT51QAULvWEMYaff/4ZDx48AABMnjwZH3/8MTiOQ1paGtatWweFQtFgbQkJCUG/fv0wfPhw9OnTB2fPnqUOSupEKpWKh3iLi4uRl5entN8t9J/s7GxwHId33nkH5ubmUCgU2LFjB3Jzc5X2WqR+ULjUk5ycHBw6dAgA0KRJE8ybNw+LFi0Sy1f89ttviI2NbZAv+LS0NMyfPx9PnjwBz/NISUnB3LlzkZWVVe+vTbQXx3Fo2rQpAKCsrAzZ2dlK+93Z2dn4/fffAQAODg5YvXo1hgwZAgCIiYnBpUuX6OFIzVG41AOh9H1iYiIAYNCgQXB2doadnR2mTJkCjuOQk5OD4ODgBmnLgQMHcPfuXQAQ56qjo6Px66+/UgcldWJnZwegopR+RkaGUn4nYwyhoaHiqH/QoEFo0aIFJkyYAH19fcjlchw5ckTcgUnUE4VLPfnzzz8hk8kglUrx5ptvitsox4wZg2bNmgEAjh8/rtSnvWcpKirC/v37wRiDhYUFvvnmG5iZmYHneezduxclJSX1+vpEuwnhwhjD06dPlfawcubMGcjlckilUgwbNgwcx6F79+5wdXUFAFy+fBnp6elKea3nYYxBLpeDMUYPYbVA4VIPSkpK8M8//wAAbG1tq1zp6uDggIEDBwIAkpKScPXq1Xr94N6+fRvR0dEAgFdffRXjx4/HK6+8AgC4efMmYmNj6+21ifazs7MTt9U/efJEKb+zuLgYV65cEX9/ly5dwHEcLCwsMGDAAADA48ePERoaWq99RyaTYcaMGZg1axZOnTolbiwg1UPhUg8ePnwoTkO1b99efLoDKg5RvvXWW9DT04NcLsdvv/1Wbx2EMYazZ8+ipKQEHMdh5MiRMDExwYgRI8BxHIqKinD+/Hl6KiO1Zm1tDX19fQDA06dPlfI7Hz58iHv37gEAvL29q1T+Hjx4MKRSKRQKBc6fP6+U13ue1NRUHD16FJs3b8aqVasabAOOtqBwUTLhvvD8/HwAFaMFPb3/P07EcRy6deuGFi1aAKg49VtfB7NkMhkuXLgAALCwsECvXr3Ew2/C2sv58+fpiYzUmqWlJYyMjAAA6enpdX5QYYwhPDwcBQUFAICePXuKIyOO49CxY0fxYe2vv/5CcXFxnV7vRe24ceOGOG396quv0s2YNUThUg/CwsLAGIOenl6VKTGBtbW1eLr/4cOHiIqKqpfRQ1pamjgl5uHhAUdHRwCAk5MT2rRpAwC4deuW0p44ie4xNzeHmZkZgIpwUcaDijBVrK+vjx49elTpP02aNEGHDh0AVFzvLSz6KxtjDOfOnQPP89DX10e/fv3q5XW0GYWLkslkMkRGRgKoGC0IX+KVcRyH119/HRKJBDKZDGfPnlV6OxhjuHXrljgq6tGjh/jkZWxsLK67ZGZm4vbt2zQ1RmrF2NhYPJWfnZ2N8vLyOv2+srIyREREAKhYr/x3/5FKpeJFYQUFBQgPD6+Xz25BQQH+/vtvAIC9vT28vLyU/hrajsJFyXJycpCQkAAAcHZ2fuZNkRzHoWvXrrCxsQEAXLp0CaWlpUpvy7Vr16BQKCCVSv9zD3nPnj0hkUggl8tx7do1pb820Q0GBgawsrICAOTl5dV5mio9PR1JSUkAgNatW8Pa2rrKn3Mchx49ekBfXx+MMXHjjLLFxcXh/v37AICuXbv+px3k5ShclCwpKQmZmZkAKhbzhfnof2vWrJk4vI+Li1P68F4ulyM0NBQA0LhxY7Rr1078M47j0KFDB/GJ88aNG7TuQmpFT09PfIAqKioS1xpr6+7du+Jou1OnTuJmgcratm0rvmZYWJjSt9MzxnD58mVxI0z//v01tpq5KtF/MSVijCE6OhplZWUAgM6dOz/37+rr66N3794AKqqtKntbZXZ2NuLj4wEALi4u4tkagb29PZycnABUhBuV0yC1wXGcuMBeWlpap80pjDFERkZCLpeD4zh06tTpmX/P2tpavO4iKSkJaWlptX7NZ+F5HlevXgUAmJqa/mfdh1QPhYuS3bx5E0DFdIGHh8dzP5Qcx6FXr14wMDAAYwwXL15UargkJiaKJ6Y7dOjwnxGUkZER2rdvD6DqVAQhNSWUgJHJZHU6pS+EC1CxluPp6fnM/qOnp4cuXboAqHgwi46OVmrfyc/Px507dwAAjo6OaNmypdJ+ty6hcFEiuVyOmJgYABWL+cLI4Hnc3d1hb28PAAgNDUVRUZFS2sEYw+3bt8URlNARK+M4Dj4+PgAqnjhv3bpFi/qkxjiOEz/DPM/X6SBlSUmJeKjXzs5O3N34LF26dIFEIoFCoUBYWFitX/NZHjx4gMePHwMAvLy8xN1wpGYoXJQoPz8fycnJACpO4r9sEdDS0lJcd3nw4IH4s8og7LgxNDREu3bt/vMEKKy7GBgYAIDSOyjRHU2bNhXXJB4/flzrh5T09HQ8evQIANCqVSs0atTomX+P4zh4enqKfx4eHq60A46MMdy8eVNcxxGqA5Cao3BRorS0NHEx383N7bmL+QKJRCJuqywqKhLPx9RVWVmZOKy3traGs7PzM/+eq6uruGOt8kiHkJqwtbUVF97rsv6RmJgobgho3779C29rbdasmXgQOS4uTqnl/oX1T319fXTu3JnCpZYoXJTo3r174hPP8+aLKxNO6xsaGoIxhr///lsp4ZKdnS1uo3RxcYGlpeUz/56lpSVcXFwAVFxsRiX4SW1YWVnB2NgYQEW41OYzLFwOJizmC+uBz2NiYiLugHzy5InSdltWPmdjY2MjXpFBao7CRUkYY4iNjQXP85BIJC9czK+sdevW4oJoRESEUspZ3L9/Xyxb4enpKU59/ZuBgYHYQXNycpQ6LUd0R+PGjcVyQk+fPq31tvbbt28DqJjKbdu27Qv7D8dx4m7M0tJSpR0ETk9PF6/KcHNzE0f2pOYoXJSEMSYu5hsbG4ulwV/G0tJS/IK/f/8+Hj58WOd2xMbGiielhTWd5xG2e5aXl+POnTu0qE9qzMTERDxImZmZWasDwWVlZeJivrW1tTjl9Twcx8Hb21s8TCnsMqure/fuidupO3bs+MxzNqR6KFyUpLy8XKzkamVl9Z9zJc8jlUrRvXt3ABUlJ5Sxa0vYDq2vr//S7dAeHh5iWZioqKg6vS7RTQYGBrC1tQVQcUpfKDpZE7m5ueKDVYsWLWBhYfHSn2nVqpX4927evAmZTFbj161MWMwXpuY6duxYp9+n6yhclCQvLw8pKSkAKvbGP2+ny78JpWD09PTAGBMPb9VWTbdDOzk5iU+dMTExda4NRXSPVCoVtyMXFhbW6iBlSkqKOJXbtm3b507lVmZjYyN+viuPOGpLqGgOVJwDq866KXk+Chclefz4sfjhdnNzq1bnELi7u4tzu2FhYXXatVXT7dBWVlZiB01OTq638v9Euzk4OACoGMHX5obI+Ph48XNfuVTRi1Q+CJyZmSluYqmt0tJSsYq4ra3tS6fmyItRuCjJvXv3xLlmDw+PGv1s5eqv9+7dq9NBtMePH4unpFu3bv3S7dBCJQEAyMrKqrcS5kR7cRwnhotcLhcPIFaXUDaJMQapVFrtzTAAxKmrsrKyOi/qZ2ZmilNzlafcSO1QuCgBYwzx8fFgjEEikaBNmzY1Gk7r6+uLp+VzcnLEaa3aSEhIELdDP+vw5LMIi/7Coiot6pOaat68uXiQ8tGjRzX6DPE8L44YzMzMxO3xLyNsWRZmCYQtxLWVlJQknpfx8vKqcskfqTkKFyWovFPMxMTkuYcWX6R79+5iOYvr16/X+qxAdHS0uB26umdtPD09xV0xwmYAQmqiWbNm4pe8cMq+uoqKisTNME2bNq1yLfjLuLi4iFO/t27dqvWaoVAySSaTidUrSN1QuChBWVmZuDfe2tpaPLdSXRzHwcvLS9wEEBoaWquzAsJBNKAi5Fq1alWtn3NxcRHL7wsH2QipCRsbG5iYmACoWJzneb7aP5ueni6e7G/VqpX4e6qjcgWK5OTkOh0EFh6sDA0NaTFfCShclCAvLw+pqakAarZTrDJ7e3uxk8TExNRqYb2kpES8qKxJkybV3g5tY2MjLl4mJiYqtZQG0Q0WFhZiJYjHjx/XaASRnJyMwsJCABWHfl9U9uXfDAwMxEX9rKysWlf3rlwyqXJ/ILVH4aIEqamp4n0obm5utTp4ZWJiIi5OPn36VAyJmsjMzBRDzsXFRTw1/TJGRkZo27YtACAjI0PcUk1IdRkbG4sXeGVkZFS7wrcwpSycLanNiKHyQeDanhOrvJnF2dn5uSWTSPVRuChBUlKSuFNMuMSoNoRLicrLy2t1N/j9+/fFUYe7u3u1FyQrzzGXlJSImxMIqS59fX1xx1heXp54ZqU6hBGDoaGhuGuyuoQpZeEgcG1P6lfehu/l5UUn85WAwqWOhHIrtd0pJhBOBAsFAK9evVqjL3ihHcIpZS8vrxq9tpeXF6RSaZWDZIRUF8dx4qVaJSUl1a6OXF5eLpZ9sbS0rNV0lJOTk1gh4NatWzUuPyOsVVa3ZBKpHgqXOhK+1IHa7xQTODk5oXnz5gAqOklNy2gIhf8MDAxeWvjv3yrfn3Hr1i2l3Y9BdIdwGFcul1e7Rl5eXp44HVXdsi//ZmlpKW5euX//fq1uwxQeqJ53/xGpOQqXOiorKxO3UdZmp1hljRo1EkccKSkpNSpiKZPJxJCzsLCo8dWsdnZ24rTG3bt3xXs1CKkOjuPg7Owsjn6TkpKqNfL+d9kXYXqrJvT19cX1ypycHNy9e7dGPy8UbQWqVqwgdUPhUkfZ2dniAriTk1O1F9GfRSKRiEUsi4qKEBUVVe2psfz8fLH8haOj40vLvvybiYmJuF705MmTGp9VIMTBwUGsCFGdXVvC4WNhGqu6ZV+eRbjUSy6XIyIiokZTytnZ2WLJJGdnZ7HWHqkbCpc6evjwobiI3qZNmzotBAr32hsYGNS4iGVqamqVsi81fQKsfD9GSUkJYmJiaFGf1Iitra14Xur+/fvVqlIslGyRSqW1no4S1gyF8zFhYWE1OmeTnJwsjp7atWtXo7qA5PkoXOqAMYa7d+/WuODei7Ru3VosYhkREVGtxUmhHULZFy8vrxp3UmHHmFCdOSwsrOaNJzqtcePG4rRwamrqS7cjKxQKcZ3Q3Ny82ncgPYujo6O4Xnn79m3x3MzLCCfzhT5MZfaVh8KljoS5Wj09vRovoj+LjY2NWEjy3r174rmVlxH299flCbB169bi/v6IiAgqv09qxMDAQFzry8rKQmZm5gv/flFRkXiey97evkZlX/7N3Ny8ynplTQqwCjXJDA0N0b59e1rMVxIKlzpQKBRiwT1zc/M67RQT6OnpoUePHgAqdtJUZ92F53nxCdDMzAxubm61em07Oztx1018fHytdt0Q3SWRSMTPXlFR0Us3pKSlpYkVwFu3bl2jsi/Peu2uXbsCAIqLi6u9XllWVoZbt24BqFpKhtQdhUsdFBQUiDvFmjdvLp5QrguO49C9e3fo6emB53n8/fffL/2ZwsJCcYdMs2bNar1jzdDQUFx3ycrKQlxcXK1+D9FdwiFIuVyOhISEF37Bx8fHi9NX7du3F6sq1wbHcejSpYt47XF11yszMjLExfxWrVrRYr4SUbjUQWpqKp4+fQqgolPV5cmrMi8vL3Hd5dq1ay9dd0lLSxMPrbVt2xampqa1fu3u3buLu25CQ0NpUZ9UG8dxVS7KE7bGP4twWFeo4O3t7V3n6ai2bduKD3ihoaEoLi5+6c/cvXtXXMz39vamk/lKROFSB3FxceIHuEOHDkqbq7W1tYWnpyeAivtZXlbrKy4uTnwC7NChQ62fAIVFfTMzMwAVwUaHKUlNtGjRQjyMGxsb+9wK2zzPiwcXzczMxNp2dWFtbS1uqklMTHzpdnph44pQ16xLly51bgP5fxQudRAZGQme5yGVSpXy5CXQ19dHr169AAC5ubkICwt77giCMSa2QxlPgC1atBAva4qKinrpoiwhlVlbW4u7thITE59bZaKgoEAc2TRv3hz29vZ1fm09PT2x3+Tn57905M3zPG7cuAEAMDU1pcV8JaNwqSW5XC4WyWvUqJFSnrwEHMehV69e4vzxhQsXnttJFAqF2A5zc/M6t8PU1FTcUPDkyRNERkbS1BipNmNjY3HdJT09/bm7HR89eiReh+zp6VmnqVyB0G+Ec2IXL1584Wc3Ly9PXMx3cHCgMvtKRuFSS+Xl5WjUqBGaNWuGVq1a1ansy7O0a9dOfJq7evXqc58A8/LyxFswW7RooZQnwAEDBkAqlUIul+P06dN1/n1EdwgHGoHnV9hmjOHWrVviORgfH586LeZX5unpKY6c/vnnnxeWMaq81b9Tp051qq5B/ovCpZZMTEywd+9ehIWFISgoSGmL+QJra2txDjg5Ofm5O7eSk5PF7Zzt27evczs4jkO3bt3Ei8ZOnz5dp9v9iG7hOA7e3t6QSqXgef65JfCvXbsGxhj09fXRpUsXpU1HWVpaiiWU7t+/L55D+zfGGK5fvy5ulunVqxdNiSkZhUsdGBgYoGnTpmjVqpXSP5gSiQQDBw4Ex3EoKSnBuXPnnvkEGB4eLp7MF3Z61ZWdnR369+8PoGLenEYvpCbatm0rLuoLC+aVlZSU4Pr16wAqDg3X9A6XF5FIJBg0aBA4jkNZWRnOnDnzzKkxnudx5coVABUPisoMOFKBwkVNcRyHV199VSxBfvr0abFEhYAxJp6DMTIyEov31ZVEIsHEiRNhbGwMhUKBPXv21PiODKK7mjVrJh5GjI2N/c/FYSkpKeL5sHbt2ol3sSiD0G+Ewq2nTp165pbkrKwshIaGAqgoOCscHibKQ+GixpycnODt7Q2gorzLv0uJ5+fnizXAHBwcan0y/984jkPXrl3xyiuvwM3NDX369KFFfVJtxsbG8PHxAVBxZbewJghUPBDduHFDLPbas2fPat+YWl2Ojo7o1q0bACA6Olosjlm5DWFhYeJ6S69evWi9pR5QuKgxAwMDDBs2DEBFkISEhFTpJLGxsWINJR8fH7EirTIYGRlhx44duHLlChYuXCjekElIdbz66qvgOA4ymQyXLl0SP7eMMZw9exaMMRgYGKB3795Kn47S09PDm2++CYlEgpKSEhw+fPg/4XLy5EnI5XJIpVK88cYbNCVWDyhc1BjHcXjjjTfEkhRHjx4VD0syxnDu3DmUlpaC4zgMGDBAaTtuhNd2cXGBra0tdTxSI8KmEKEI6rlz58Qp3czMTFy+fBkA0LJlyxpdx12T1x80aJC4c/Lo0aPiphegYkpMWEds3rw5unXrRp/xekDhouZcXFzQu3dvABUVmIWnwOLiYpw4cQJAxQ4Z4UmREHXg6OgoTunevHkTsbGx4tkT4eR8v379anWtcXXY29tj1KhRAIAHDx5g//79YIyBMYaQkBCxntjrr7+ulJqA5L8oXNScnp4e/P39oa+vj/Lycmzbtg0lJSW4evWqeACsR48eNb7WmJD6ZGBggFGjRoHjOBQWFuLgwYMoLS3F3r17oVAoYGhoiLfffrveHog4jsOUKVNgYWEBxhi2bNmCpKQkZGdnY8uWLeB5HiYmJvDz86OHsnpC4aLmOI5D//79xQXK8+fP4/vvv8f69etRVlYGPT09+Pn5KX1RlJC64DgOw4YNEw8XBwYGYuvWrTh//jyAiku56nM6iuM4tGvXDuPHjwdQcebl448/xqxZs8SzN4MHD6YtyPWIwkUDmJqaYt68eTAxMYFMJsP8+fNx5swZABUL+YMHD6YOQtROixYt4O/vD47j8PjxY3z++ecoLS2Fvr4+AgIClFLy5UWkUinmz58vFrM8d+4cgoODwRhD06ZNsXjxYrrSuB5RuGgAYYFy5syZ0NfXh0wmA2MMtra2WLFihXhgjRB1wnEcZs6cKV7ipVAoIJFIMH78eIwcObLeH4g4joODgwN++OGHKvfFODo6YsuWLUqtZE7+i+ZSNIS+vj4WL14MZ2dnHDlyBObm5pg+fTqVrSBqi+M4NG3aFPv27cO6desQGxuLvn37YsaMGTAyMmqwNvj4+ODMmTO4cOECSktL8eqrr8LZ2Zn6TT2jcKkFxhhiY2NhaGjY4K/drl07uLu7g+M4SCQS8ZRxfUtKShK3lhLCGENMTEy1+4Cfnx8UCgWkUini4+PruXXP5uTkBKDi9snaXOGdnJxcb7vbtBGFSy107doVf/31FyIiIlTdlAbD8zxeeeUVVTeDqIkuXbrg77//fm5hSm1EfaBmOEZ1PWpE1/9z0VQCoT5AfaA6KFwIIYQoHe0WI4QQonQULoQQQpSOFvQ1ROXZS5rzJbqK+oHmoJGLhoiMjISenp5O7c4h5N8iIyMhkUioH2gAChdCCCFKR+FCCCFE6ShcCCGEKB2FCyGEEKWjcCGEEKJ0FC6EEEKUjsKFEEKI0lG4EEIIUToKF0IIIUpH4UIIIUTpKFwIIYQoHYULIYQQpaNwIYQQonQULoQQQpSOwoUQQojSUbgQQghROgoXQgghSkfhQgghROkoXAghhCgdhQshhBClo3AhhBCidBQuhBBClI7ChRBCiNJRuBBCCFE6ChdCCCFKR+FCCCFE6ShcNABjDDk5OQCAnJwcMMZU3CJCGh71A81C4aLGcnNzsXHjRri5uWHAgAHgeR4DBgyAm5sbNm7ciNzcXFU3kZB69+9+AID6gQbgGMW/Wjp9+jRGjx6N4uJiAKjylMZxHADAxMQER44cweDBg1XSRkLqG/UDzUXhooZOnz6NoUOHgjEGnuef+/ckEgk4jkNISAh1LKJ1qB9oNgoXNZObmwsHBweUlJS8sEMJJBIJjI2NkZKSAgsLi/pvICENgPqB5qM1FzXz888/o7i4uFodCgB4nkdxcTH27t1bzy0jpOFQP9B8NHJRI4wxuLm5ISkpqUY7YTiOg4uLCxISEsR5aEI0FfUD7UDhokYyMzNha2tbp5+3trZWYosIaXjUD7QDTYupkcLCwjr9fEFBgZJaQojqCGdZaov6gXrQU3UDyP8zMzOr08+vX78evXr1go+PD5ydnWlqgGiEkpISxMfHIzo6GrGxsbh9+3adfp+5ubmSWkbqgqbF1Igw15yYmFjjn7WyskL//v0RExMDhUIBKysrdO7cGT4+PujSpQtat24NqVRaD60mpGby8/MRGxuLmJgYxMTEiGsrjRs3hoeHBzw8PDBt2jQ8fPiwRr+X1lzUC41c1AjHcRg+fDg2bNhQ459bsmQJAgICUFRUhKioKISHhyM0NBTr1q1DeXk5zMzM0KlTJ/j4+MDHxwft2rWDvr5+/fxDCKkkKytLDJKYmBg8evQIANCkSRO4u7tj8ODB8PDwQLNmzcRQ+PTTTzF79uwaLegzxhAQEEDBoiZo5KIm5HI5tm7dit27dyMuLg4KhaJaHetl+/vLy8tx+/ZthIWFITw8HBERESgqKoKhoSE6dOiAzp07o0uXLvD29oaxsXE9/MuILmGM4cmTJ4iOjkZMTAxiY2Px9OlTAICDg4M4MvHw8ICNjc1zf09Nz7kAgJ6eHsLCwtChQwel/FtI3VC4qIFHjx5h/vz5iIuLwyeffIKmTZti+PDh1T6ZfPLkSQwaNKhar6VQKBAXFyeObMLDw5GTkwOpVApPT09xZNOpUyc0btxYWf9EoqUYY3jw4EGVkUlubi44joOzszM8PDzg6ekJd3d3NGrUqEa/uyYn9AFg9OjRkEql8PX1xdChQ2kEo2IULip24sQJrFixAtbW1lizZg08PT0BVL+m0tGjR6sdLM/CGENSUhLCwsLE/z158gQA0Lp1azFsOnfujCZNmtT6dYh2UCgUSExMFBffY2NjUVRUBD09Pbi5uYmjkjZt2sDExKTOr1eTftC3b1/s378fv//+O7y9vTF9+nR6QFIhChcVKSoqwooVKxASEoLhw4fj888/h6mpaZW/k5ubi71792LTpk1VFvldXV0REBCAiRMnKr3zMMbw+PHjKmFz//59AECLFi3EsPHx8YGDgwM9HWq5srIyJCQkiNNcd+/eRVlZGYyMjNCmTRsxTNzc3GBgYFAvbahpP7h58ya2bNkCAJg+fTq8vb3rpV3kxShcVOD27dtYsGABcnJysHjxYgwZMuSFf58xhuzsbBQUFMDc3BxWVlYN+qWemZmJ8PBwMWzi4+PBGEOTJk2qhI2rq6s4RUE0U1FREeLi4sQprnv37kGhUMDMzAzu7u7w9PSEh4cHnJ2dG3z3YU36QV5eHrZt24bIyEgMHToU48ePpw0sDYzCpQHxPI+ffvoJW7duhYeHB1avXg0HBwdVN6vGCgoKEBERIa7Z3L59GwqFAo0bNxan0Hx8fODh4UHbn9Vcbm4uYmNjxZHJgwcPwBiDlZVVlcV3R0dHjRulMsZw8uRJBAcHw9HRETNnzoS9vb2qm6UzKFwaSHp6OhYuXIiwsDC89957mDp1KvT0tGMneGlpqbj9OSwsDFFRUSgtLYWxsTE6duwojmy8vLxgZGSk6ubqLMYYMjIyqiy+P378GADQtGnTKmFiZ2encWHyPMnJydi4cSOysrIwefJk9O3bV2v+beqMwqUBXLhwAUuWLIGhoSFWrlyJLl26qLpJ9UomkyEmJgZhYWEIDQ1FREQECgoKoK+vj/bt24ujm44dO9a5KgF5PsYYUlNTERMTIy7AZ2ZmAgBatmwpTnO5u7vDyspKxa2tX2VlZdizZw/OnTuHHj164MMPP/zPGidRLgqXelRWVoZ169bhwIED6Nu3L5YuXaqTd00oFAokJCSIZ21CQ0ORlZUFiUQCd3d38axNp06dtP5Lrj4pFArcv39fnOaKjY1Ffn4+JBIJXF1dxVGJu7u7zob6tWvXsHPnThgbG2PmzJlo06aNqpuktShc6sm9e/cwf/58PHr0CJ999hnGjBlDQ/H/YYzh4cOH4sgmPDwcKSkpACp2AAkjmy5duqBp06Yqbq36kslkuHfvnjjFFRcXh5KSEujr66NNmzbiyKR169Y0HVlJZmYmNm3ahPj4eLz99tt46623aG2wHlC4KBljDAcOHMC6devQokULrFmzBq1atVJ1s9TekydPqoxshC2nzZs3F4PGx8cHLVu21NmQLi0tRXx8vDjNlZCQAJlMBmNjY7i7u4sjE1dXV9oZ9RIKhQJHjx7F4cOH0aZNGwQEBLywYgCpOQoXJcrNzcWSJUtw8eJFjB07Fp9++ikMDQ1V3SyNlJOTI+5ICwsLQ2xsLHieh7W1dZXtz25ublr71FlQUPCfAo88z6NRo0ZikHh6eqJly5a0BbyW4uLisGnTJpSUlOCjjz5C9+7dVd0krUHhoiQ3btzAwoULIZPJsGzZMvTp00fVTdIqhYWFiIqKEs/a3Lp1CzKZDObm5lUKcnp6emrsU3t2dnaVnVxCVWAbGxtx4d3T0xP29vY6O3qrD0VFRdi1axeuXr2K/v37Y9KkSfRQqAQULnUkFJz86aef0KVLF6xYsYLKpDSAsrIysSBnWFgYIiIiUFJSAiMjI3To0EEMG29vb7VcbxAKPAojk+joaLHAY/PmzatsC67LrYykehhjuHDhAn766SdYW1tj5syZcHZ2VnWzNBqFSx08evQICxYsQGxsLGbMmIFJkybR9ISKKBQKxMbGVilbk5eXB6lUinbt2lUpyFnTAorKIGxiqDwyycnJAcdxcHJyqrKTSxd3FKqL1NRUbNy4ESkpKfD19cWQIUNolFhLFC61FBISghUrVsDS0hJr1qxBu3btVN0kUgnP80hKShJ3o4WFheHp06fgOA5t2rQRqwj4+PjUy0KuQqFAUlKSGCSxsbEoLCyEVCr9T4FHOm+hXmQyGfbt24eQkBB07NgR06ZNowKYtUDhUkOVC04OHToUixYtoi8HDcAYQ0pKihg0oaGh4ppGy5YtxRs7O3fujObNm9f4abW8vBx3794Vz5jEx8ejrKwMhoaG/ynwSPP5miEqKgpbt24FAMyYMYPuiakhCpcauHPnDubPn4+cnBwsWrQIQ4cOVXWTSB1kZGSIU2jh4eG4e/cuGGNo2rRplbM2Li4u/wmb4uLiKgUeExISoFAoYGpqWmWKy9XVVWt3s+mCvLw8bNmyBTdv3sSwYcMwbtw4jd0w0tAoXKqB53ns2bMHW7Zsgbu7O1avXg1HR0dVN4soWX5+PiIiIsSRTXR0NBQKBSwtLdGuXTvY2dnBwMAA2dnZePjwIRhjsLCwECsFe3h4oEWLFjRHr2WEAphBQUFo0aIFZs2ahWbNmqm6WWqPwuUlMjIysHDhQoSGhmLKlCmYNm2a1hScJM+XmZmJiIgInD17FteuXcP9+/eRn58PjuNgbm6Odu3aoXfv3ujfvz+8vLxoqksHJCcnY8OGDcjOzsaUKVPQp08fepB4AQqXF7h06RK+/PJL6OnpYdWqVejatauqm0TqgXBBmjDFFR0dLRZ4dHR0FM+YtGrVSqwkIGx/LiwshIGBAby8vMQqAt7e3rQOp6VKS0vx008/4cKFC1QA8yUoXJ6hrKwM69evx/79+9GnTx8sW7aMtodqEZ7n8eDBA7G4Y3R0tFjg0cXFpcqaibm5+XN/j0KhwN27d6tsf87OzoZUKoW7u3uV7c+WlpYN+C8k9e3q1avYuXMnTE1NERAQQAUwn4HC5V8SExMxb948PHz4EJ999hneeecdGvpqOLlc/p8Cj8XFxdDX14ebm5s4MmnTpg2MjY1r/TqMMdy/f79K2Aj3pbRq1apK2Ro7Oztl/fOIimRkZGDTpk1ISEgQC2DSObf/R+HyP4wxHDp0CN9++y0cHBywZs0auLm5qbpZpBZKS0tx9+5dMUzi4+PFAo9t2rQRF+BbtWpV7zt/0tLSqoRNUlISAMDBwaFK2NBGAM2kUChw5MgRHDlyBG3btsUnn3xCBTD/h8IFFQUnly1bhvPnz+Pdd9/FnDlzaIFWgxQWFoplVGJjY5GYmAiFQgFzc/MqZVScnJxUvi04KytL3JEWFhaGuLg48DwPGxsb8ZyNUJCTnoI1R2xsLDZt2oTS0lJ8/PHH6Natm6qbpHI6Hy5hYWH4/PPPUVZWhmXLlqFv376qbhJ5iZycnP8UeGSMwdraWqwU7O7uDgcHB7UfDRQUFCAyMlI83Hnr1i3I5XKYm5uLQdOlSxd4eHjQLkU1V1RUhJ07d+LatWtUABM6HC5yuRzbt2/HDz/8AB8fH6xcuZIKTqohxhjS09PFxfeYmBikpaUBAOzt7f9T4FHdw+RlSktLxYKcoaGhiIqKEgtyent7i6ObDh06qGVBTl3HGMP58+fx008/wdbWFrNmzULLli1V3SyV0MlwSUlJwYIFCxAdHY0ZM2Zg8uTJNAWhJhhjePToUZWRSXZ2NjiOQ8uWLauMTHRhB59cLkdMTIxYRSAsLAz5+fnQ09ODl5eXWEWgY8eOL9zZRhpWamoqNmzYgNTUVPj5+eGNN97Q+AefmtK5cDl58iS+/vprWFhYYM2aNfDy8lJ1k3SaQqFAcnJylTARCjy6urqKi+9t27al8wSo2EZ979498cbOsLAwZGRkiAU5K6/bWFtbq7q5Ok3XC2DqTLgUFRVh1apV+P3336ngpArJZDIkJCQgOjpa3MlVWloKAwOD/xR4pGmflxNGepVrpAkFOZ2dnatcEW1vb6/i1uqmyMhIbN26FRKJBNOnT9eZApg6ES53797FnDlzkJmZiUWLFmHYsGGqbpLOiYqKwqFDh5CQkAC5XA4TE5P/3PtOC9bK8fTpU3EKLSwsDAkJCQCApk2bokePHlixYoXOTdGoWm5uLrZs2YJbt25hxIgR8PX11fr3QCfCRSaT4cmTJ2LhQdLw8vLy8PTpUzRq1Ajm5uYwMTHR+s6lLhQKBUpLS1FSUgKe52njioowxpCWlgaZTKYTi/w6ES6EEEIaVoPOQ+h6jqnDkzq9B/QeqBq9B6rXEO9Bg09yX7lyBbGxsWrxAWsojDF0795dbRbywsPDkZSUpHPvQYcOHdSmwOClS5cQExOjc+/BK6+8Am9vb1U3BQAQERGB5ORkVTejwbVv3x6tW7eu99dp8HD566+/YGNjAxcXl4Z+6QYlk8mQkZGBZs2aiSev1SVcIiMjYWFhAQcHB2RkZMDS0lLrF9Pv3LmD+Ph4tQmXK1euwMbGBq6urqpuSr3ieR5Pnz6FnZ0dwsLCEBUVpTbhUrkfaDOZTIa8vDzY2NggOjoad+/e1c5w4TgOnTt3RufOnRv6pRsEYwxlZWX48ccfcfLkSSxfvhzu7u4oKChQddOq8PDwQHl5OU6dOoV+/fph1KhRWh0wJSUlKCwsVHUzRBzHiVuEtRFjDPn5+QgKCsKvv/6Kr7/+Gh4eHsjPz1d100Qcx4m7FbWVTCbD4cOHce3aNXz00UdwdnZGcXFxg7w2HUtXIsYYiouLsXHjRvzwww9ISUnBihUrkJ2dreqm/UdGRgZ27dqFp0+f4vDhwzh16hQUCoWqm0W0AM/zuHXrFqZPn47t27fj8ePH+O6779TuAUvbyeVy/Prrr/jtt9+QlpaG77//Hnl5eQ32+hQuSsIYQ15eHlauXIkDBw5AoVDAzMwMvr6+ankq18bGBsOHD4eBgYF4kvjixYsUMKTWGGMoKSlBUFAQpk2bhoiICPA8DysrKwwePFinizg2NLlcjt9//x3Hjh2DQqGAkZERRowY0aDfRRQuSsAYw9OnT7F48WKcPHkSPM/DwsICX375JUaOHKmW000SiQSDBg3CmDFjoKenh7KyMuzZswfXrl0Dz/Oqbh7RMEIf+PLLL/Htt98iJycHEokE3bp1w86dOzFp0iQ6Y9ZAFAoFzpw5g8OHD0Mul8PQ0BD+/v7o3bt3g24gUb9vPQ3DGENycjKWLFmC27dvAwCaNWuGL7/8Et26dVPrgph6enoYPnw4SktLcezYMRQXF+P777+HkZEROnXqpFM7mUjt8TyPqKgofP3114iLiwMAmJmZYeLEiZgwYQLMzMzos9RAeJ7HpUuX8Msvv0Amk8HAwADjxo1Dv379Gvy7SH2/+TQAYwy3bt3CnDlzxGBxcXHB2rVr0b17d7UOFoG+vj5Gjx6NwYMHQyKRoKCgANu3b0d0dLTOnwUgLyeTyXD8+HHMnDlTDBZnZ2d8++23+Pjjj2Fubk7B0kB4nsf169exd+9elJWVQU9PD2+//TYGDRqkkkvy1P/bT03xPI8rV65g7ty54l75Dh06YN26dWjXrp1GdSgDAwOMHz8effr0AcdxyMnJwdatW5GQkEABQ56JMYaioiJs3boVX331FbKysiCRSNC7d29s374dr776qspv/dQljDHcvn0bP/zwA4qLiyGVSjF8+HAMHTpUZdPyFC61ICyWLV68GOnp6eA4Dr1798Y333wDJycnjQoWoGJLppGRESZNmoTu3buD4zikp6dj8+bNePDgAQUMqYIxhszMTCxZsgQ//PADSktLYWhoiIkTJ2Lt2rVwdHTUuD6gyRhjSEhIwI4dO5Cfnw+O4zBw4ECMHj0a+vr6KmsXhUsNMMZQXl6OoKAgrFq1Cvn5+ZBKpRgxYgS++uorjb4JkeM4mJiY4IMPPkDHjh0BAI8fP8bmzZuRlpZGAUMAVPSB+/fv49NPP8Uff/wBhUIBCwsLLFq0CDNnzqRpsAYmXLmwbds2ZGVlgeM49OrVC+PGjVNpsAAULtXGGENpaSm2bt2KrVu3ineQTJw4EQsWLNCKTsVxHBo1aoSpU6fC09MTAHD//n1s2rQJGRkZFDA6TlhjDAgIQHh4OBhjcHR0xLfffou33nqLdoM1MMYYMjIyxLNEANCxY0dMmjQJRkZGKv8+onCpBsYYCgsL8c033yAoKAgymQzGxsYICAjA1KlTYWxsrPI3Ulk4joOlpSVmzJiBVq1aAQASEhKwZcsWZGdnU8DoKJ7ncfHiRcyaNQuJiYkAAC8vL2zevBk9evTQiM0r2kQ4V7djxw7x/XB3d8eHH36oNrvz6BPxEsKb+PXXX4sHkho3boyFCxdi7NixKh961geO42Bra4uAgAA4OTkBAKKjo7F9+3bk5+dTwOgYuVyOY8eO4fPPP8fTp0/BcRz69OmDjRs3ws3NTS2+yHSJUAlk9+7diI6OBgC0bNkSU6dOhaWlpdq8HxQuL8AYQ3Z2NpYuXYozZ86A53nY2Nhg2bJlKt2F0RA4joO9vT0CAgLQvHlzABW3SW7fvh2FhYUUMDqivLwcgYGBWLFihbjG+Oabb2LVqlWws7NTmy8yXSKs+4aGhoIxBjs7O0ybNk3t3g8Kl+cQ5jMXL16MixcvgjGGZs2aYdWqVejdu7dOTANwHIcWLVogICAAdnZ2YIwhLCwMu3btQlFREQWMFhPWGHfs2IGNGzeipKQE+vr6mDJlChYuXIhGjRqp1ReZrpDJZDhy5AguXLgAxhgsLS0xdepUtdylqv3fkLXAGEN6ejoWL16Mq1evAgAcHR2xevVq+Pj4qN2bWJ84joOrqysCAgJgY2MDxhiuXr2KH3/8ESUlJRQwWkiYdlm/fj12796N8vJyGBsbY+bMmZg+fTpdUa0iCoUCf/zxB06cOAGe52Fqaor3338f7u7uavl+ULj8i1AjaeHChbhx4waAihPHa9asQfv27dXyTaxvHMehTZs2+OSTT2BlZQXGGC5fvoyff/4ZpaWlFDBaRNi8smrVKvzyyy+Qy+UwNzfHwoULMWHCBNoRpiJCWZdDhw6J9cImTpyIzp07q+13EoVLJZULUIaHhwOoKOeyevVqtX06aCgcx8HT0xPTpk2DhYUFGGM4f/48goKCUF5eTgGjBYTNK8uXLxc3r1haWmLZsmVaf9+POhOmowMDA8WyLmPGjMGrr76q1tPz6tuyBiassXzxxRcICwsDALi6umLNmjVo3bq1TgeLgOM4eHt7izWjeJ7Hn3/+if3796O8vFzVzSN1wBhDbm4ulixZIlb2trGxwYoVK1RWm4pUvC/R0dH4/vvvUVRUBKlUimHDhmHIkCFq/55QuOD/y1l88cUXCA0NBVARLKtXr0arVq0oWCoRbhL96KOPYGpqCoVCgZCQEBw+fBgymUzVzSO1IOyKXLx4Mc6ePQvGGJo2bYo1a9bozOYVdcQYQ2JiIrZv3468vDxwHId+/fph9OjRGjGK1PlPDWMMOTk5WLZsGa5fvw6gYo1l5cqVFCzPIdzT8cEHH8DY2BgKhQK//fYbfv31VwoYDSMEyxdffCHuirS3txcre9PnXzUYY0hNTcXWrVuRkZEBjuPQo0cP+Pr6asy6l06Hi3DP99dff42///4bANCiRQusXLmSpsJeQiKR4JVXXsHkyZNhZGQEuVyOo0eP4sSJE5DL5apuHqmGysFy6dIlsZzLN998o9YLxdpOmKLfunUrUlNTAVRUXJ8yZYpGVQPR2XARSoavXbtW3DPevHlzrFy5Em3bttWYN1CVpFIp+vTpgwkTJsDQ0BAymQwHDx7EmTNn6LpkNSeM2L/88ssqwbJ27Vp4e3vT519FhPdl27ZtYlmXNm3a4KOPPtK4+oU6Gy6lpaXYuHEj/vjjD/GU61dffQVPT0+NegNVTSqVYsCAARg/fjz09fXF08MXL16kgFFTwq6wpUuXilNhQrDo6nZ7dcAYQ0FBAXbt2oWYmBgAFWVdpk2bBisrK417X3QyXMrLy7Fz504cPXoUPM/DysoKS5cuRceOHTXuDVQHUqkUr7/+OsaMGQM9PT2UlZVhz549uH79OnieV3XzSCXCF9hXX32Fc+fOiSP21atXU7CoEGMMJSUl+PHHHxEZGQmg4rr0GTNmoGnTphr5vuhcuMhkMgQGBiI4OBgKhQKNGjXC4sWL0a1bN418A9WFnp4eRowYgZEjR0IqlaK4uBi7du1CVFQUnYFRE8LJ+5UrV+LUqVPirrDVq1fTVJgKMcbEB7KrV6+CMQZbW1tMnz4dLVq00Nj3RafCRaFQ4NixY9i1axdkMhlMTEzw2WefoU+fPrTdUgn09fUxevRoDB48GBzHoaCgANu3b0d8fDwFjIoJT8bffvstTpw4AcYYmjRpgpUrV6JTp04a+wWm6YQLCPft24fLly9XqRem6RWndeYbled5nDt3Dhs3bkRZWRkMDQ0xffp0DB06lIJFiQwMDDB+/Hj07t0bHMchOzsbW7ZsoeuSVay8vBxbtmzB4cOHwfM8rK2t8dVXX9GIXcVkMhkOHTqEP//8EzzPw9zcHB999JFWrP3qxLcqYww3btzAqlWrUFhYCKlUikmTJuGdd95R+1OumobjOBgZGWHSpEnw8fEBAKSlpWHLli14+vQpBYwKlJeX4/vvv0dQUJB4H9GSJUvQs2dPjf8C02QymQy//vorQkJCoFAoYGpqig8++EBrpii1PlwYY4iNjcXy5cuRk5MDiUSCMWPGYPLkyVp50Zc64DgOZmZm4hMYACQnJ2Pbtm3Izc2lgGlAcrkc+/btw+7duyGXy2FqaoqFCxeiX79+NGJXIblcjt9++w3Hjx+HQqGAsbExJk2ahK5du2rN+6Id/4rnYIzh4cOHWLJkCR4/fgyO4zB48GDMmDEDhoaGqm6eVuM4DhYWFpg2bRqcnZ0BVNxmSXfBNByFQoHjx49j8+bNYtn8zz77DEOGDNGaLzBNJJfLceLECRw9elSscOzv76/2hShrSnv+Jf8i1Atbvnw5EhISAADdu3fH3LlzYWpqqhXDTnXHcRzs7OzE7ZQAEBoair1796KsrIwCph7xPI+zZ89i7dq14kVf06dPx+jRo2kqWIXkcjn++OMPHDp0CDKZDIaGhvD19UXfvn21KlgALQ0X4U6KNWvWiKXzPTw88MUXX6jVHdO6gOM4tGzZEtOnTxfvgrlw4YJ4LwVRPsYYrl+/jq+//hoFBQXQ09PD5MmT4efnpxEFD7WVXC7HqVOncODAAchkMhgYGGDs2LEYOHCgVga+VoZLeXk5Nm/ejPPnzwOouEVy6dKlaNasGQWLCnAcB3d3d3z44YcwNTUFz/M4ceKEuJBJlIcxhjt37uCLL75AVlaWuMb48ccfa0zBQ20kBItwPYW+vj7Gjh2LwYMHa2WwAFoYLnK5HHv37hVP31tbW2PJkiUav2dc0wml+idPngxDQ0PI5XIcPHgQFy9epFP8SsIYw4MHD7Bo0SJxjfH111/H7NmzaY1RheRyOU6ePFklWN599128/vrrWj2S1Kpw4XkeISEh+OGHH8SdMfPnz6cKr2pCIpHgtddew9ixY8UyMT///DPCwsJo/aWOhEq6ixcvxr179wAAPXv2xMKFC2mNUYWExfsDBw5UCZYhQ4ZodbAAWhQujDH8888/WL9+PUpLS2FgYIDp06ejf//+1LHUiFQqxRtvvIERI0ZAKpWiqKgIO3fuRHR0NAVMLQlXRyxbtgwREREAAC8vLyxdupTWGFVIOMciLN4bGBhg3LhxOhEsgJaEC2MMMTExWLFiBfLy8iCVSuHr64sxY8Zo7XymJtPX18fbb7+Nvn37guM45ObmYuvWrUhMTKSAqYWSkhJ88803uHjxIgDAyckJK1asoDVGFRFKuhw6dEi8QM/Q0BDjx4/X+qmwyjQ+XBhjSElJwbJly5CWlgaO4/DGG2/ggw8+0Jk3URMZGBhgwoQJ4m2H6enp2Lx5M1JSUihgakCo8H38+HGxXtjXX38NV1dXChYVEIpQ7tu3D7///jvkcjmMjIzg7++PwYMH69R3kkaHi3CxzvLly3H37l0AQLdu3TBnzhyNurFNF3EcBxMTE3zwwQfo0KEDACAlJQWbNm2iMjHVpFAosH//fvz8889ihe8lS5bQ1REqIlSd3rNnD06dOgWFQgETExNMmTIFAwYM0LlZFI0NF+GNXLt2LcLCwgAAbdu2xRdffAELCwvqXBqA4zg0atQI06ZNQ9u2bQEASUlJ2Lx5MzIzMylgXoDneZw8eRKbNm0ST9/PmzdPLBhKGpZwT87OnTtx4cIFsQjlhx9+iNdee03rDkhWh8b+i2UyGbZt24YzZ86IFx4tXboU9vb21Lk0CMdxsLKywowZM+Di4gIAiIuLw9atW5GTk0MB8wzC5pXVq1ejuLgYenp6+PjjjzFixAid/BJTNcaYWP37+vXrYIzBwsICU6dORffu3XX2PdHIf7VcLkdQUBAOHDgAnudhaWmJL774Am3atKFg0UAcx6Fp06YICAhAixYtAAC3b9/Gtm3bkJeXRwFTiXBIcsmSJcjJyRE3r0ycOFGn5vPVBWMMaWlp+O6778SL8WxtbTFz5kx07txZZ4MF0MBw4Xkev//+O3bt2gW5XA4TExPMnTsXXbt2pWDRYBzHwcHBAQEBAbC3twcAREVFYfv27cjPz6eAQdVDkpU3r8yYMYNO36sAYwzJyclYv369uObr4OCA2bNnw8PDQ+e/jzQqXBhjuHjx4n/OsgwePFinnxC0BcdxcHJywsyZM9G0aVMwxhAeHo4dO3agoKBApwOGMYb09PT/HJL8/PPPYWJiouLW6R7GGKKjo7F+/Xo8fPgQANCqVSt8+umntFPvfzTmG1m48GvFihUoKCiAVCrFhAkT6MIvLcNxHFxdXREQEIAmTZqAMYbQ0FDs3LkThYWFOhkwwiHJ5cuX/+eQJG1eaXg8z+P69evYuHEj0tPTwXEc2rdvj9mzZ6N58+b0fvyPRoQLYwy3b9/G0qVLxWJ8b775Jt577z268EsLcRyH1q1bIyAgADY2NmKVX10NmH8fknR2dqZDkiqiUChw/vx57NixA3l5eeA4Dt27d8cnn3wCGxsbej8qUftwYYwhPj4eixcvFueZBw4ciJkzZ8LIyEjVzSP1hOM4tG3bFjNnzhQD5tq1a9ixY4dOBUx5eTm2bdsmHpK0s7OjQ5IqUl5ejuPHj2PPnj0oLi6GVCrFwIED8dFHH6FRo0b0fvyLWocLYwyJiYlYuHChOK/Zq1cvLFiwAGZmZvRmajmhVH/lgLl+/Tq2b9+uE2swQoXvwMBAKBQKNG7cGEuWLNGaO9Y1BWMMpaWl2L9/Pw4dOiQWoBw1ahQmTJhAB7afQ23DhTGGpKQkfP7550hKSgJQcfr+yy+/pHlmHfKsgLlx4wa2bdum1bvIFAoFjhw5gu3bt0Mmk8HY2Bjz58/Ha6+9Rp/9BiRcPLh7926cPHkSCoUCRkZG8PX1xejRo2FgYEDvx3OoZbgIU2Hz5s0Tryju3Lkzli1bRvOaOkgImFmzZomL/GFhYdi8ebNWHrTkeR5//PEH1q1bh5KSEhgYGCAgIADDhg2jXZENSLgqfdOmTbhy5UqVU/e6VICyttTuk8oYw82bNzFv3jwkJiYCqAiWr7/+GnZ2dhQsOqryGoydnR0YY4iMjMSGDRuQnp6uNQHD8zwuXLiAFStWoLCwEHp6enj//fcxbtw4+jJrQMKZonXr1uHmzZtgjMHGxgaffPIJXnnlFQr5alCr/0IKhQIXLlzA3LlzxTWWrl27YsWKFWjatCkFi47jOA5t2rTB7NmzxYOW0dHR+Pbbb5GcnKzxAcMYw19//YUlS5aIV0eMHz8e77//Ph2SbEA8z+PmzZv45ptvxCl5R0dHfPrpp+jQoQMFSzWpxX8lxhhKSkrw888/44svvkBGRgY4jkOfPn2wcuVKGrEQEcdxaNWqFebMmQMnJycAFcUu16xZg+vXr0OhUKi2gbXEGMOVK1ewaNEiZGdnQyKRYPTo0QgICKBdkQ2EMQaZTIYzZ85g48aN4veQp6cn5s6dSzv0akjl42ye55GcnIwtW7bg8uXLUCgU0NPTw8iRIzFz5kyYm5vTG0qq4DgOLVu2xNy5c7Fjxw7cuXNHnBsfOHAgRo4cqVE3MPI8j0uXLuHLL78Uz3GNHDkSc+fOhbGxsaqbpxOEqsYHDhzAhQsXIJfLIZFI0LNnT0yYMIG2GteCysKF53lkZGTgt99+w4EDB5CZmQkAMDMzwwcffIB3330XhoaG9IaSZ+I4DnZ2dpg9ezaCg4Nx6dIllJWVISQkBOHh4RgwYAA6d+4MKysrtZ5SUigUOHPmDL766ivk5uaKwSKUdaHPf/3jeR6JiYnYu3cv7t69C8YYDA0NMWLECAwfPpy+h2pJJeFSVFSEjRs34vTp0+LFUJXn0318fKikC3kp4T6Y999/H61bt8bBgweRlZWFtLQ0BAYG4vDhw2jcuDHat28vlvNXJ3K5HIcPH8a6detQWFgIqVSKN998E/PmzYOpqSl9oTWA8vJynDx5Er/++isKCgoAANbW1pgwYQK6du1K30N1oJJw4TgOERERePLkCQDA0tISb731FsaPHw8rKyvqVKTaOI6DgYEB+vfvDw8PD/z222+4du0aCgoKUFJSgpKSEjg7O6vlImx6ejp++OEHcVfY2LFjMXPmTBqxNKDi4mKcO3cOBQUF4vrKpEmT4OjoSO9BHakkXExMTDBq1ChkZ2djwIABGDVqFFq2bAmO4+gNJbXCcRyaNWuGDz74AMOHD0dkZCRiY2ORkZEBJycntfxcNWvWDAsWLMDy5csxZswYcVeYOrZVW1lYWGDs2LH48ccfMWDAAAwdOpRO3CtJg4eLcPLeyckJ8+bNg6WlJQoLCxEdHd3QTWkwKSkpaNSokaqbUcWjR4+0+txE69at4ebmJu4eu3nzJkxNTVXcqv8n9ANhKtjBwQFxcXGqbla9Urd+wBhDSkoKWrRoAT8/P9jZ2SE1NVXVzapXT58+hZmZWYO8VoN/u7Rr1058qtQVjDF4e3uruhmiVq1aIS4uTtzDrwsYY2jbtq2qmyHy8vJCREQEYmJiAEAspa/NeJ5Hx44dVd0MkdAPkpOTAUB8L7RZQ/YDjjXgyTNNP+RWV+ow1Kb3gN4DVaP3QPUa4j1o0HAhhBCiG9RvCw0hhBCNR+FCCCFE6XQiXGQyGR49eoTy8nJVN0Vn5ebmIi4uDjzPq7opOquoqAhPnz6l90BFGGN4/PgxHjx4oOqmNAidCJfk5GRMmzYNvXv3xu+//67zi3kNSSaTYc+ePZgyZQoOHjyI4uJiVTdJZ12+fBlDhgzBG2+8gdu3b6u6OTolJycHK1aswKxZs/DXX3/pxHeQzizoFxUVYfXq1fjtt98wZMgQLF68WK3OPWij1NRUrF+/Hg8fPsSECRMwbNgwtdgppMsePnyIuXPnIjo6GrNmzcKUKVPUsnqBNomMjMTWrVshkUgwY8YMtG/fXtVNahA6Ey6CP/74A1999RUsLCywevVqnXmjGxJjDOfPn8fu3bthY2ODTz/9FM7OzqpuFvkfuVyOTZs2Yffu3ejWrRvWrFmDJk2aqLpZWkcmkyE4OBgnT55Ep06dMG3aNLU6RFrfdC5cgIon6gULFuDOnTuYPn06Pb0pUVFREXbs2IG///4bAwYMwJQpU+g+EjV17do1zJs3DzKZDCtXrkTfvn1V3SStkZqaig0bNiA1NRX+/v54/fXXdW7UrpPhAlQ8ve3YsQO7d+9G586dsWrVKnp6q6O4uDh89913KCoqwtSpU9GzZ09VN4m8RE5ODhYtWoQLFy7A19cXc+fOhaGhoaqbpbGEUftPP/0EW1tbzJo1Cy1btlR1s1RCZ8NFEBYWhs8//xxlZWVYtmwZPb3VAs/zOHz4MA4cOCDWyrK1tVV1s0g1Mcbwyy+/YM2aNXBycsL69evh6uqq6mZpHGHUfv36dQwYMAATJ07U6aDW+XABKrbJLlu2DOfPn8c777yDzz77TKc/FDWRmZmJDRs2IDY2Fu+88w7efvttugNDQ929exdz5szBo0ePsGDBArz77rs6N5VTW7Gxsdi0aRNKS0vx8ccfo1u3bqpukspRuPwPYwyHDx/GN998AwcHB6xZswZubm6qbpZau3r1KrZt2wZjY2PMmjULHh4eqm4SqaPS0lKsXbsWv/zyC/r374+vv/4aFhYWqm6W2lIoFDh8+DCOHj2Ktm3bIiAgANbW1qpullqgcPmXxMREzJ8/Hw8ePMCcOXPo6e0ZSktL8dNPP+HPP/9Ejx49MHXq1AYr400axrlz57Bo0SIYGRlh7dq16Nq1q6qbpHYyMjKwadMmJCQkYMyYMXjzzTdpY1AlFC7PUFZWhu+++w6//PILevfujeXLl9PT2/8kJydj/fr1yMzMxHvvvYf+/ftT+GqpJ0+eYP78+QgNDcVHH32E6dOna/UdQDXxzz//YNeuXTA1NUVAQADatGmj6iapHQqXF7h06RK+/PJL6OnpYeXKlTo9j8oYQ0hICPbu3QtHR0d8+umnaN68uaqbReqZQqHA7t27sXnzZrRr1w7ffvstHBwcVN0slRFG7RcuXMArr7yCDz74gA5jPweFy0tkZGRg0aJFuHHjBiZPnqyTT295eXnYvHkzIiIiMGzYMPj7+0NfX1/VzSIN6ObNm/jss8+Qm5uLpUuXYujQoapuUoNLTk7Ghg0bkJ2djffeew+9e/emUfsLULhUA8/z2LNnD7Zs2QJ3d3esXr0ajo6Oqm5Wg4iKisKmTZvA8zwCAgLQqVMnVTeJqEhBQQGWL1+OEydOYNSoUTpTQkkYtQcHB6NFixaYNWsWmjVrpupmqT0Klxq4c+cOFixYgOzsbCxatEirn97kcjmCg4Nx/PhxeHt7IyAggNadCBhj+O2337B8+XLY2Nhg3bp1aNeunaqbVW9yc3OxdetW3Lx5E8OHD8e4ceN0buaitihcaqioqAgrV67EiRMnMHToUCxatEjrnt4eP36M9evX48GDB/D398fw4cNp+E+qePjwIebMmYPY2FjMnj0bkydP1rqdUpGRkdi2bRsAYMaMGejQoYOKW6RZKFxqKSQkBCtWrIClpSVWr14NLy8vVTepzhhjuHDhAnbv3g1LS0t8+umndFKbPJdMJsPGjRvxww8/4JVXXsHq1au1ojKDTCbDvn37EBISgo4dO2LatGlo3LixqpulcShc6iAlJQXz589HbGwspk+frtFPb0VFRdi5cyf++usv9OvXD++//z4VnCTVcvXqVcybNw8KhQKrVq1C7969Vd2kWktNTcXGjRuRkpICPz8/vPHGGzRqryUKlzqSy+XYtm0bfvzxR/j4+GDlypUaVwCTCk6SusrOzsbChQtx6dIl+Pn5aVwJJaHg5J49e2BtbY1Zs2bByclJ1c3SaBQuSnLjxg0sXLgQMpkMy5YtQ58+fVTdpJfieR5HjhzBgQMH0Lp1a8yaNUvjgpGoD8YYgoOD8c0338DZ2Rnr1q3TiGlVYdR+7do19O/fH5MmTdKoYFRXFC5KJJwBuHDhAt59913MmTNHbT+klQtOjhkzBmPGjKGCk0Qp4uPjMWfOHKSmpuLzzz/HmDFj1HZqKS4uDps2bUJJSQk++ugjdO/eXdVN0hoULkrGGMPBgwexbt06ODo6Ys2aNWjVqpWqm1XFtWvXsG3bNhgaGmL27NlUcJIoXWlpKVavXo0DBw5g4MCB+Oqrr9RqUVyhUODo0aM4fPgw2rZti08++QQ2NjaqbpZWoXCpJ/fu3cP8+fPx6NEjzJkzB++8847Kn97Kysrw008/4cyZM+jevTumTZtGBSdJvfrzzz+xePFimJiYYO3atejSpYuqmyQWnLx7965YcJJG7cpH4VKPysrKsH79euzfvx99+/bF0qVLVXYQ8f79+1i/fj3S09Px3nvvYcCAASoPO6Ibnjx5gnnz5iE8PFwsgKmqL/OrV69i586dMDExwcyZM6ngZD2icGkAFy9exJIlS6Cvr49Vq1Y16NMbYwwnT57E3r170bx5c3z66ac6XXiQqIZCocCuXbuwdetWeHl54dtvv23QwqfCqP38+fPo0aMHPvzwQ607/KxuKFwaSHp6OhYuXIiwsDBMmTIF06ZNq/cyEvn5+di8eTPCw8MxdOhQTJgwgQpOEpWKiorCZ599hry8PCxbtgxDhgyp99dMTk7Gxo0bkZWVhSlTpqBPnz40am8AFC4NqHIBTA8PD6xevbreRhE3b97Exo0bwfM8PvnkE3Tu3LleXoeQmiooKMDSpUtx8uRJvPnmm+KajLIJo/bg4GA4Ojpi5syZsLe3V/rrkGejcFGBO3fuYP78+cjJyVF6AUy5XI59+/bh2LFj6NChAwICAmBpaam030+IMjDGcPz4cSxfvhxNmjTBunXr4OnpqbTfn5eXh61btyIqKgrDhg3DuHHjaNTewChcVKSoqAgrVqxASEgIhg0bhoULFz53DpgxhqysLBQWFsLMzAzW1tbPHNanpaVh/fr1uH//Pvz8/DBixAga/hO1dv/+fXz22WeIj4/H7NmzMWnSpOeWUKpuP7h58ya2bNkCAJg+fTq8vb3r859AnocRlTpx4gTr3r07Gzp0KLtz506VP8vJyWEbNmxgrq6uDID4P1dXV7ZhwwaWk5PDGGOM53l2/vx5Nm7cODZt2jR27949FfxLCKmd8vJy9s0337C2bduy9957j2VkZFT58+r2g/Lycvbzzz+zMWPGsBUrVrDc3FwV/GuIgEYuauDRo0dYsGABYmNjMWPGDEyaNAl//vknRo8ejeLiYgAVT20C4WnNxMQEwcHBSEpKwpUrV6jgJNFo//zzD+bPnw+e57Fq1Sq89tprOH36dLX6wa5du3Dz5k08fPgQfn5+GDJkCI3aVYzCRU3I5XJs3boVP/30E+zt7XH69GkAFZsAnkfoPH369MHy5cvRq1evBmkrIfUlKysLCxcuxOXLl9G9e3fs3bsXQPX6wciRI7F+/Xo4Ozs3SFvJi1G4qJmzZ89i8ODBL+xM/2ZiYoLU1FS6KZJoBcYYdu7ciWnTpqG6X08cx8HY2Jj6gRrRzMtHtFh0dHSNggUASkpKxCc8QjQdx3EoKyurdrAAFYFE/UC90MhFjTDG4ObmhqSkpBp1LI7j4OLigoSEBJpnJhqP+oF2oHBRI5mZmXW6JjYzMxPW1tZKbBEhDY/6gXagaTE1UlhYWKefLygoUFJLCFEd6gfagcJFjdS1/L25ubmSWkKI6lA/0A4ULmrE2toarq6uNZ4v5jgOrq6usLKyqqeWEdJwqB9oBwoXNcJxHD755JNa/WxAQAAtYhKtQP1AO9CCvprJzc2Fg4MDSkpKqrUlWSKRwNjYGCkpKbS/n2gN6geaj0YuasbCwgJHjhwBx3HPLeAnkEgk4DgOR48epQ5FtAr1A81H4aKGBg8ejJCQEBgbG4PjuP8M84X/z9jYGCdPnsSgQYNU1FJC6g/1A81G4aKmBg8ejJSUFGzYsAEuLi5V/szFxQUbNmxAamoqdSii1agfaC5ac9EAjDFkZ2ejoKAA5ubmsLKyokVLonOoH2gWChdCCCFKR9NihBBClI7ChRBCiNJRuBBCCFE6ChdCCCFKR+FCCCFE6ShcCCGEKB2FCyGEEKWjcCGEEKJ0FC6EEEKUjsKFEEKI0lG4EEIIUToKF0IIIUpH4UIIIUTpKFwIIYQo3f8Ba+XfBguSIxsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 500x400 with 7 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot(beta=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64d2573b",
   "metadata": {},
   "source": [
    "### Fix the first layer activation to be linear function, and the second layer to be sine functions (caveat: this is quite sensitive to hypreparams)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e2e78752",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best value at boundary.\n",
      "r2 is 0.9982730203177389\n",
      "Best value at boundary.\n",
      "r2 is 0.9978885507047336\n",
      "Best value at boundary.\n",
      "r2 is 0.9967572093934152\n",
      "Best value at boundary.\n",
      "r2 is 0.9977250982977139\n",
      "r2 is 0.9737003463841577\n",
      "r2 is 0.9826072748927763\n"
     ]
    }
   ],
   "source": [
    "for i in range(2):\n",
    "    for j in range(2):\n",
    "        model.fix_symbolic(0,i,j,'x')\n",
    "        \n",
    "for i in range(2):\n",
    "    model.fix_symbolic(1,i,0,'sin')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fae3f32",
   "metadata": {},
   "source": [
    "### After setting all to be symbolic, we further train the model (affine parameters are still trainable). The model can now reach machine precision!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "308b72af",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pde loss: 4.38e-15 | bc loss: 5.92e-16 | l2: 4.99e-16 : 100%|█| 20/20 [00:09<00:\n"
     ]
    }
   ],
   "source": [
    "train()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35985ae9",
   "metadata": {},
   "source": [
    "### Print out the symbolic formula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f0ec310e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.5*sin(-3.14159*x_1 + 3.14159*x_2 + 7.85398) + 0.5*sin(3.14159*x_1 + 3.14159*x_2 + 4.71239)]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula, var = model.symbolic_formula(floating_digit=5)\n",
    "formula"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
